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ABSTRACT 

The particle mass used in cosmology N-body simulations is close to 10 10 Mq, which is about 
10 65 times larger than the GeV scale expected in particle physics. However, self-gravity in- 
teracting particle systems made up of different particles mass have different statistical and 
dynamical properties. Here we demonstrate that, due to this particle mass difference, the 
nowaday cosmology N-body simulations can have introduced an excessive core collapse pro- 
cess, especially for the low mass halos at high redshift. Such dynamical effect introduces an 
excessive cuspy center for these small halos, and it implies a possible connection to the so 
called "small scale crisis" for CDM models. Our results show that there exist a physical limit 
in cosmological simulations, and we provide a simple suggestion based on it to relieve those 
effects from the bias. 
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. 1 INTRODUCTION 
t-H ■ 

■ With the rapid development of computer science, the N-body cos- 
' mology simulation has become an important method for studying 
, dark matter particle systems. Such numerical experiment method 
' showed the Cold Dark Matter(CDM) universe with a dark energy 

\l ' parameter A can have nice agreement with observations on large 
£-f) . scale topics (e.g. Springel et al. 2006, Coles 2005). But the nu- 
f^-) merical predictions on small scale topics depart from the obser- 
f^*) , vations: High-resolution rotation curves of low surface brightness 
. galaxies show the halo density profiles have flat cores (e.g. Burk- 
s' ' ert et al. 1995, de Blok et al. 2002;2005,Gentile et al. 2005), yet 
• i-H , the simulation results tell us they should have cuspy centers (e.g. 

■ Navarro et al. 1997;2004); Simulation results also predict about 10 
' to 100 times more subgalaxies round our Galaxy(e.g. Klypin et al. 

. 5r ! 1999, Moore etal. 1999a, Springel etal. 2008). This is the so called 
"small scale crisis" and caused people's suspicion on CDM models. 

Many different explanations were carried out to displace 
the traditional CDM models, such as the Warm Dark Mat- 
ter(WDM)(e.g. Colombi, 1996), the Self Interacting Dark Matter 
(SIDM)(e.g. Spergel et al.2000, Dave et al. 2001), the Modified 
Newtonian Dynamics (MOND)(e.g. M. Milgrom, 1983) models 
and etc.; yet the new models have also caused new disputes for 
themselves(e.g. Yoshida et al.2003, Markevitch et al. 2004, Zhao et 
al. 2006, Kuzio de Naray et al. 2010). Anyway, the N-body simula- 
tions have brought strong support and also challenge for the CDM 
models. Notice that a lot of projects based on simulation dark mat- 
ter halo properties are now on going, for instance, the Weak Inter- 
action Massive Particles (WIMP) direct detections(e.g. H. Wang, 
1998) and the various discussions about galaxy formation et al., so 
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one can see that people do need reliable simulation results that can 
greatly help us understand the mysterious dark matter. 



2 PHYSICAL BIAS 

Here we notice that there exist a physical bias in nowaday simula- 
tions: the particle mass. Due to the technical limitation, the particle 
numbers used in cosmological simulations are limited. To obtain 
the mean density of the universe , people have to set a huge mass 
for each particle in simulations. 

With the improvement of computer science, the particle num- 
bers used in simulations have increased from 10 5 to about 10 9 
within the last decade, and the particle mass used on small scale 
topic has decreased from 10 10 Mq to lO 3 A/0 (e.g. Springel et al. 
2008). It is acceptable to set the particle mass as lO n M when 
studying the evolution of large scale structure, for we can explain 
each particle as one galaxy. But for the small scale topics, such 
as the dark matter halo property, the galaxy formation, the galaxy 
merging process, the first star formation and etc., the simulated par- 
ticle mass is still about a factor of / ~ 10 65 times larger than the 
expected GeV candidates in particle physics(e.g. Gaitskell et al. 
2004). At the same time, the dark matter number density also has 
the same factor smaller than expected (that means using only about 
10 1 ~ 10 s particles to simulate one 1O 13 M dark matter halo). If 
we don't think the mass of one dark matter particle can be heavier 
than our human body, we should consider that a physical bias exist 
in simulations: 

m sim — > rriDM x /, 

Tlsim — > TlDM I f (1) 
Do these two kinds of self-gravitating particle systems with such a 
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bias have the same statistical and dynamical properties? If not, one 
should be cautious when applying the simulation results. 

2.1 Long Term Relaxation 

It is not a simple question to describe the difference from the bias 
in general. But for a virialized dark matter halo, the numerous theo- 
retical and numerical studies on globular clusters can give us much 
help. One point is the long term relaxation effect. 

From the view point of one particle, when it is flying in a stable 
and spherical halo potential, theoretically its energy Ei and angular 
moment Li should be conserved as constants. But for a virialized 
dark matter halo, the potential is contributed by many moving par- 
ticles, that means the potential will no longer be ideally spherical 
and stable. In this case, both Et and Li can be changed. Intuitively, 
when the halos include fewer particles, such effect will be more 
serious. Note that such long term relaxation effect we attention is 
caused by particle density field fluctuations on large scale, but not 
by collisions of a few close particles. 

Binny & Tremaine (1987;2008) show the relaxation time scale 
treiax caused by long term particle density field fluctuations (here 
the softening parameter does not affect the results) should be: 

t retax ~ (0.lN/lnN)t cross (2) 

where N = M/rrii is the particle number of the halo, t cros3 = 
R/v is the crossing time scale and v is the virial velocity (v 2 ~ 
GM/R) . Analytical and simulated results (see Huang S. et al. 
1993, Diemand et al. 2004) give the similar formula. For one halo 
with given M, t r£ i ax oc Inf / f. The Eq.lO shows us that the bias 
greatly shortens the relaxation time. 

If we define the mean free path as L 3 = t re i a x/v, then we 
can follow the SIDM models (e.g. Spergel et al.2000, Dave et al. 
2001) defining the "scattering cross section" as a = l/(Lp). In 
the central region of a typical simulation halo, the scattering cross 
section is about a S i m ~ 9 x 10 -26 cm 2 /GeV (see Xiao et al. 
2004), that is approximately the value expected in SIDM models 
((Tsim — 0-Ictsidm)- In contrast, for the GeV CDM particles 
ccdm — 10~ (Tsidm — 0. 

Now we find the difference: the bias has bring an excessive 
scattering cross section for the CDM models. The value of a s im 
cannot be neglected for CDM models, but not big enough for SIDM 
models. Will it affect the dynamical properties of the halos? 

2.2 Core Collapse 

The excessive scattering cross section means particles in simula- 
tions will have an excessive way to exchange their energy and an- 
gular momenta. Then the simulation halos are possible to follow 
the evaporation effect appearing in globular clusters: Once a par- 
ticle exchanges its energy and gets Ei > 0, it can fly away and 
never come back. In a virialized system, the mean particle energy 
< Ei >= —GM/R < 0. That means the evaporating particles al- 
ways bring out energy, and the left particle system becomes tighter 
and tighter. Such process appear more serious at the central part of 
the halo, and the result is to introduce a dynamical core collapse of 
the system. 

Such evaporation and core collapse processes have been well 
studied in galactic dynamics on the topics of stellar clusters. Since 
the dark matter halos in nowaday simulations are similar to the 
globular clusters: both are virialized systems and consisted of pure 
gravitational interacting particles, and even have the similar parti- 
cle numbers (about 10 1 to 10 s ); we can use the same method to 



estimate their core collapse time scales. Following the way analyz- 
ing stellar clusters (see Spitzer, L. 1969, Giersz et al. 1994, Binney 
et al. 1987) we get the core collapse time scale of a virialized dark 
matter halo in simulations: 

t ~t ( M \h m r x r Th n\ 

tcc-M^^) l 3xlO 8 M0 i 4ofcp C j "^ 

Here M is the halo mass, and m is the particle mass, is the half 
mass radius, t u = 1.37 x 10 10 yr is the Hubble time in ACDM 

models. The parameter S = , A 9 . , -r-J— ~ 10° is not sensi- 

tive to the particle numbers of the halo Nhalo, and we have sug- 
gested t cc to be about 16 times of the half mass relaxation time t r h 
(see Takahashi 1995). 

Before discussing in detail, we should emphasize the effect of 
the softening parameter e introduced in simulations. Softening is 
a numerical trick introduced in N-body simulations to prevent nu- 
merical divergences when two particles become very close (and the 
force goes to infinity), the method is to modify each particle gravi- 
tational potential, such as the form $ = 1 . The introduc- 

Vr 2 +E 2 

tion of e can effectively affect the short term "two-body relaxation" 
process. However, 

(1) The softening parameter e is unable to make the halos 
avoid such core collapses. Because the gravitation is a long term in- 
teraction, the relaxation process discussed above is mainly caused 
by the long tern particle encounters. The introduction of e has no 
business with these long term process. In fact, the time scale deriva- 
tion of eq.l(2} and eq.10 in Galatic Dynamics is based on the dis- 
cussion of the density field fluctuations in distance and e will not 
change it. 

(2) One other point is that e prevents the hard binaries for- 
mation. The hard binaries release energy and drive a reexpansion 
of the core after the core collapse in a globular cluster (e.g. Cohn 
et al. 1989), yet the softening parameter e makes such processes 
impossible for dark matter halos in a cosmological simulation. 



3 A PHYSICAL LIMITATION 

Equation I0 can tell us many secrets. For a given dark matter 
halo (with setting value of M and r^), we find the core col- 
lapse time scale is proportional to the particle number of the halo: 
tc, oc N oc 1/m. For the GeV CDM particle halo, t cc > t u 
and the core collapse will never happen within one Hubble time. 
But for one Galaxy dark matter halo in simulations, if we use less 
than N ~ 10 12 Af o /3 x 10 8 Af Q ~ 3000 particles to progress the 
simulations, the bias of particle mass will bring an excessive core 
collapse within one Hubble time. Our result shows a limitation of 
the particle numbers ~ 10 3 that should be used when studying the 
Galactic scale topic in simulation. 

How about the CDM halos in cosmological simulations? In a 
given simulation (with setting particle mass m), Equation l(3j shows 

1 1 

us tcc oc M 2 r£ . This means the larger dark matter halos will 
have core collapse later, or to say, smaller halos at high redshift 
will be more "dangerous". Since the ACDM models show us a 
hierarchical structure formation scenario, we expect to avoid such 
core collapse process in the whole cosmological simulation, if we 
ensure all the smallest halos at the beginning follows the limit: 

tcc ^ at u (4) 

The at u (a ^ 1) is the mean time scale of these smallest halos 
existing in the universe before merging. If we set a = 1, then we 
ensure the core collapse process caused by the relaxation effect will 
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not happen in the smallest halos (so for all the larger halos within 
the whole hubble time). 

The parameter of the smallest halos at high redshift are de- 
cided by the initial conditions of the simulations. In nowaday 
cosmological simulations, people apply the linear theory and use 
Fourier power spectrum P(k) to describe the initial fluctuations 
S(x), and to generate the initial conditions(e.g. Seljak U. et al. 
1996, Springel et al. 2008). But due to numerical limitation, the 
simulation initial conditions can only represent part of P(k) in a 
limited range [kmin, k max ], where k m in is decided by the simula- 
tion box size, and k m ax figures the smallest halo properties at the 
beginning. 

In hierarchical structure formation scenario, the smallest halo 
formed by the collapsing of the dark matter within one shortest 
wavelet A = 2n/k max . So we estimate the mass of it as: M ~ 
p^(±) 3 ~ Ppr 2 -. In halo models (see Gunn et al. 1972), the 



spherical collapse halos have the mean density of about 178p, then 
their characteristic radius ro follow 178pro ~ p(^)' i , if we set 
r/j ~ O.laro, (for NFW density profile 2 ^ a ^ 3 ), then we get: 



r h ~ 0.1^=4 _ 

VV78 2 

find the limit: (here po 



5.58 X 10 



— 2„ 



. Combining eq.((4j and eq.([3]l, we 

= 1.4 x W^Mq/Mpc 3 ) 



m < to* 



68M Q ( 



2tt 



kmax X WkpC po 



(5) 



The parameter m* is a physical limit. Wether the Virialized dark 
matter halos show us the dynamical difference within one Hubble 
time, depends on the relation between the used particle mass to and 
the limitation m« . When we use particles with mass m ^ m* in the 
simulation, we can ensure all halos avoid the core collapse process 
caused by the unexpected relaxation effect. But if one use heaver 
particles in simulations (m > ra„), the hugely magnified scatter- 
ing cross section can introduce an excessive core collapse for all 
the smallest halos, and the cores do not re-expand like a globular 
cluster due to the softening parameter. The dynamical difference 
between the "ant particle" (m ~ GeV) system and "elephant par- 
ticle" (m ~ 10 65 Mq) halo will be serious in this case. 



3.1 An Excessive Core Collapse 

The limitation given by eq.l[5) seems easy to fulfil, for it is also 
much larger than the GeV scale. In a AC DM universe with p ~ 
po, if we generate the initial condition with kmax = 27r/(73 x 
lOfcpc), then we calculate the m, should be 2.7 x 10 7 Mq. This is 
the case in Navarro's simulation (1997). Unfortunately, the particle 
mass m used in this work is 6.8 x 10 9 M©, that's much larger than 
the limitation, to > m* implies that an excessive core collapse 
have been introduced in that simulation, and we find similar results 
appear in many popular works. 

We can understand the limitation to* in another way. In a 
simulation with given particle number N and box size L, one can 
rewrite eq.lO as k max fc 3 oc m _1 p2 ; where to = L 3 p/N is 
the particle mass used. Then we can express the limitation with the 
Nyquist frequency k ny — nNs /L as: 



Lbi 

where b is a factor: 



= k 



ny i , 
7T&3 



(6) 



6 = 0.44(1 + z)ln(0.1N halo ) a - * (-£-) 3 . 

Po 

For one halo with 100 particles at red shift z ~ 5, the value of 



b is close to 1, and its cube root implies b is not a sensitive fac- 
tor here. The Nyquist frequency is corresponding to two particles 
within one wavelength, so a small halo of fc* scale includes about 



(tt x 2) 3 



250 particles. Now we can understand the physical 



limitation: we should ensure the smallest halos (corresponding to 
wavelength 2n /k m ax) include at least 248 particles. All the halos 
with particle numbers less than this limitation will be dangerous in 
the simulation processes. 

However, people set the k ma x to be the Nyquist frequency 
in Navarro's simulation(1997), so their smallest halos include only 
about 2 3 = 8 particles and the limitation in eq.l[6) is violated. Al- 
though we are able to use much more particles (about 10 3 times 
more) in the recent simulations, when generating the initial con- 
ditions to represent P(k) in [k m i n ,k max ], people are still artifi- 
cially setting k ma x to be the technical limitation k ny , rather than 
the physical limitation fc» . Therefore, in these simulations k ma x = 
k ny = irb 3 fc* and is always larger than fc* , no matter how large the 
particle number N is, and the smallest halos are always containing 
about 8 particles, the limitation in eq.([6j is always violated. Hence 
we find these simulations have definitely introduced an excessive 
core collapse for their high redshift small halos. 



3.2 Fossils in the Simulation ? 

Though we have demonstrated the unavoidable existence of the ex- 
cessive core collapse process for the small halos qualitatively, it is 
not easy to make clear how will it tamper with the simulation results 
quantitatively. However, the dynamical effect of the high redshift 
small halos is apparent: the excessive core collapse (maybe more 
than one times for the same halo) can make the halos to be more 
concentrated than they should be, and makes their density profiles 
to be much more cuspy at central regions. 

In a hierarchial structure formation scenario, these high red- 
shift small halos will soon bring themselves into complex merger 
process: they merge into each other or are devoured by huge halos. 
A lot of theoretical and numerical studies have been carried out on 
the merger process and shown us the qualitative properties: 

• Major Merger process: When two halos with similar mass 
merge together, theoretical study (e.g. Dehnen, W. 2005) and nu- 
merical experiments (e.g. Taffoni G. et al. 2003, Michael et al. 
2004, Ileana et al. 2008) showed us the two halos syncretize to- 
gether, the center of one halo sinks rapidly to the center of the 
other halo. But the central density sloop information can be well re- 
tained: a major merger of two-cored halos yields a one-cored halo; 
yet mergers between a cuspy halo and a core/cuspy halo, the inner 
density sloop of the remnant will be closer to that of the steeper one 
of the initial systems. 

• Minor Merger process: For the merger process between a 
large halo and a satellite (M s ^ O.lMh), Taffoni G. et al.'s semi- 
analytical and N-body method study (2003) showed us the fate of 
the satellite halo is determined by its orbit and concentration prop- 
erty: low concentration satellite below O.lMh is disrupted by tides 
quickly; yet high concentration one can survive with a low mass 
center and become a new substructure of the large halo. 

Since the dark matter halos are assembled step by step from 
the high redshift small halos in hierarchical CDM halo models, one 
can image how the excessive core collapse affects the simulation 
results: 

(1) It make the earlier small halos to be too concentrated and 
leave them a much too steep density sloop at the central region. 

(2) Such too cuspy character can be retained during all of the 
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major merging process later on, including the simulation final prod- 
ucts at z — 0. 

(3) It also makes some small halos which should be disrupted 
when swallowed by huge halos become too concentrated, and leave 
an unwanted substructure after the minor merger. 

So our discussion imply a possible connection with the un- 
expect cole collapse and both the 'cuspy problem' and 'substruc- 
ture problem' : The too cuspy feature of halos and the numerous 
substructures maybe are both the fossils of the unwanted core col- 
lapses? 

Without quantitative numerical experiments, we can only get 
such qualitative suspects here. It is possible that some other physi- 
cal mechanisms bring about the unconformity between simulation 
halo properties and observations. Yet since our discussions have 
shown the excessive core collapse can bring unwanted fossils, when 
we are trying to discuss these possible mechanisms, we should 
avoid these unwanted fossils to get reliable simulation result. 

3.3 Can we Avoid it? 

It is the small halos corresponding to the wavelength between fc» 
and k ny who introduced these fossils. The reason is that the tradi- 
tional method can not give enough particle numbers to model the 
halos on this scale. Hence we suggest, in a simulation with given TV 
and L, we should ensure the physical limitation k max ^ fc, when 
setting the initial conditions, but not use the technical limitation 

kmax — kny- 

For one excessive subhalo relate to the fossil within one high 
redship small halo, yet the too cuspy density profile of a huge halo 
corresponding to all fossils within each small halos of its merger 
tree, we can expect when people use fc, < kmax < k ny in simula- 
tion, they can see the subhalo number decrease serious but the too 
cuspy density profile will not change unless they use kmax < k, . 

Some numerical experiments have already shown us interest- 
ing results: Moore et al.( 1999b) have compared results with limited 
k' max , which is less than k ny (though still larger than and the 
normal case {kmax = k ny ). Their results show the huge halo still 
contain too cuspy density profile (see figure3), but the number of 
substructures dropped seriously (see figure4). We can now expect 
the too cuspy density profile can also be changed when they use 
k 1 <r k 

"•max ■ 

If our suspects is correct, one can also make a simple "predic- 
tion": Since the dangerous small halos correspond to wavelength 
of fc, < k < k ny , the traditional method simulation (setting 
kmax ~ kny) with higher resolution (k ny will be larger, corre- 
sponding to larger number of smaller halos), will show us a much 
more larger number of smaller subhalos in the simulation result. 
Just as shown in Via Lactea simulation (e.g. Kuhlen M. et al. 2008). 
And we predict similar result can still be seen if Kuhlen M. use 
more higher resolution. 

In case that when k ma x and L have already been given by the 
physical object of the simulation, we can calculate the minimum 
particle number that is needed from eq.([6jl: 

TV > TV, = k max L s b (7) 

The traditional method has set k ny = kmax-, then the particle num- 
ber used is N ny = k^axL 3 /7r 3 , it is always about ir 3 b ~ 31 times 
fewer than what we need. Here we suggest to use TV ^ TV, parti- 
cles in the simulations and avoid the unwanted core collapse. The 
estimation of TV, upside is just a rough calculation, but we do need 
it, for modeling a 10 60 particle halo with only about 8 simulation 
particles is unacceptable in physics. 



When kmax and L are given, one can also examine the effect 
from the bias with similar analysis above: we can expect subha- 
los number decrease when TV > N ny is used; and the too cuspy 
property of large halo can be changed when TV > TV* . 

Another point is that TV, in eq.lJTJ can help us avoid the exces- 
sive core collapse, but not every thing from this bias. Some other 
topics, for instance, the translation of angular momentum, should 
be studied further in details. 



4 SUMMERY 

In summery, we discussed a long term relaxation effect which is 
amplified hugely by the bias of particle mass/number in cosmolog- 
ical simulations. With the mature theories people used in studying 
globular clusters, we find such relaxation can not be neglect: 

• We find a physical limit exist: we should use at least about 
250 particles to model the smallest halo. If not, the relaxation pro- 
cess can bring an excessive cole collapse within one Hubble time, 
especially for the small halos at high redshift. 

• Our discussion imply such unwanted core collapse process 
can leave fossils in the final halos. 

• Unfortunately, people used to set kmax = kny in simula- 
tions. It means the smallest halo include only about 8 particles and 
the physical limit are always broken. 

• We give a simple suggestion to avoid such effect. By setting 
kmax ^ fe* but not kny, we can get more reliable result. 

The dynamical properties of CDM particle systems people 
collect from N-body simulations have already been basic of many 
popular astrophysical topics. Hence we suggest people attach im- 
portance to such a dynamical effect, abating the unwanted fossils 
and get more reliable results. 
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